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ABSTRACT 


Soe SSS 


Nonlinear mathematical programming techniques and their 
application to the solution of selected optimal control 


problems are reviewed, implemented and evaluated. 


Pertinent techniques for unconstrained minimization 
problems are discussed and applied to a set of test 
problems. The same procedure is followed for a ‘set “of 
constrained problems. Performance data is given for all 


test problems. 


Two means of transforming the optimal control problem 
into a mathematical programming problem are outlined. The 
mathematical programming techniques are applied to the 
solution of two optimal control problems, the Earth-Mars 
rocket problem and the minimum energy regulator subject to 
linear constraints. The relative performances of the 


techniques are then evaluated. 
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CHAPTER I 


In the last two decades there has been considerable 
research and progress in the fields of optimal control and 
mathematical programming. Many techniques based upon the 
Calculus Ofe (Vartations [17], sDynamic wProqgranuningoel 44, 
Pontryagin's Maximum Principle [32], and their extensions 
have been developed to cope with the numerical solution 
problem. But it is only in the last several years that 
direct mathematical programming techniques have been applied 


to the numerical solution of optimal control problems. 


The aim of this thesis is to review, implement and 
evaluate standard nonlinear mathematical programming 
techniques and then compare their relative performance in 
the solution of certain mathematical programming and optimal 


control problems. 


42. Hustorical Backqround 


eee a a ae Se i we SS ee 


The mathematical programming (MP) problem involves the 
extremization of an objective function subject to inequality 
or equality constraints. The classical technique used to 
solve the equality constrained optimization problem is 
Lagrange'ts Method of Undetermined Multipliers. Even though 
some success has been achieved in adapting this method to 
problems involving inequality constraints, the trend has 


been to evolve computational methods of attack specialized 
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to unconstrained problems, notably the various gradient 


techniques. 


In 1941 Courant [6] proposed another method of handling 
LHe TCONEtralnts of functional minimization problems, now 
known as the method of penalty functions. This work was 
basically ignored until the early 60's when it received 
renewed interest. At this time gradient-type computational 
metheds were develceped, and, more recently, improved methods 
such as conjugate-gradient and second-variational methods 
have cone into being. "Exact" methods of handling 
constraints including such techniques as Rosen's Gradient 
Projection 33%] and Bryson and Denham's State Space 
Constraint Technigue [8] have also been developed. The 
extension of penalty methods includes Fiacco and McCormick's 
Sequential Unconstrained Minimization Technique (SUMT) 


[9,10,11]. 


In the design of optimal control systems the classical 
method is the Calculus of Variations, an alternate view of 
which is Pontryagin's Maximum Principle. Another method 
which has found wide conceptual application is Bellman's 
Dynamic Programming {4]. Lately, however, the techniques of 
Mathematical Programming have been applied to the solution 
of optimal control problems, of particular value being the 
paper by Lasdon, Warren and Rice [26] which extended the 
application of nonlinear programming methods to a wider 


class of nonlinear control systems as compared to previous 
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authors. Recently there have been several papers applying 
Fletcher-Reeves conjugate-gradient method [14] and Davidon- 
Fletcher-Powell's variable HeETich .algorathaie[ 7,494 sco 
cont Tole problens’ [ 25,27, 31,36,43:).. iImiaddition Balakeishnan 
has developed another technique called the Epsilon Method 
[2,3], which lends itself very readily to solution by 
mathematical programming techniques as will be shown in this 


work. 


Chapter II reviews the more pertinent techniques for 
solving unconstrained minimization problems, restricting 
itself to those technigues employing first order gradient 
information but with a brief mention of a second order 
technique. The algorithms for these methods are given and 
experimental results are presented at the end of the 


chapter. 


Chapter III reviews the better known techniques for the 
optimization of constrained minimization problems with 
emphasis on the SUMT and Gradient Projection (GP) 
algorithms. The algorithms for these two methods are given 


with experimental results. 


Chapter IV reviews some of the methods for solving the 
optimal control problem. + discusses briefly the methods 
of Calculus of Variations, the Minimum Principle, and 


Dynamic Programming. The discrete formulation of the 
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continuous optimal control problem and its solution by means 
of mathematical programming is discussed. Finally there is 
a discussion of Balakrishnan's Epsilon Technique and its 


solution using a Rayliegh-Ritz expansion. 


Chapter V applies the methods implemented in Chapters 
Ii and III to the solution of certain control problems by 


means of the techniques discussed in Chapter IV. 


In Chapter VI there is a summary of the results of 
previous chapters. Conclusions are drawn and 


recommendations for further study are proposed. 
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2. Numerical Optimization Strategies 


Se eS SS 


In this chapter the methods of solving unconstrained 
minimization problems are reviewed. Emphasis is placed on 
descent methods which make use of first-order gradient 
information. Direct search techniques such as Rosenbrock's 


method are ignored. 


The unconstrained minimization problem consists of 
finding the vector x*=(x*,x#,...,xX*) which minimizes a 
scalar funceion. f(x}. This unconstrained minimization 
problem is a special case of the Mathematical Programming 
problem (see Chapter MIII). The special feature of it is 
that the solution need not satisfy any constraints or 


conditions. 


Section 2.7 of the chapter deals with techniques for 
minimization along a line. These are methods for solving 
optimization problems in one variable. Such methods are 
necessary for all of the search techniques dealt with in 


this chapter. 


Zoutendijk [45] gave the name “methods of feasible 
directicns" to descent methods that generate a sequence of 
feasible directions. This category of methods, as defined, 


includes most of the search procedures in common use. 
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The general feasinia cicccr on method is defined as 
follows: 
For small “, by a first order Taylor series expansion, given 
a direction of search d: 


f (x+«d) =f (x) + «< VE (x) ,d> 


And for a direction of descent 4d 
£(x+«d)<£ (x) 
iS Gsadtisfied: d is: *usable" at x if 
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Additionally, for constrained problems a must be 
sufficiently small that no constraints are violated. 
(tSetect nvtiall trial point «x* =xre xX; 
(2) Generate a feasible direction a“ corresponding to x‘; 
(3) Find the value of %,>0 minimizing £ (x*+«,a‘) and obtain 
the next trial point as 
xkti=xk + aa‘; 
(4) Set k=k+1 and repeat from step (2) as required. 
Cost_Factors: 
This procedure solves the original vector minimization 
problem by generating a sequence of optimizations whose 
efficiency is dependent upon the algorithms used to 


generate A and d. In determining the cost of a 


particular algorithm one must consider the following: 


(1)ease of programming 
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(2) storage requirements 
(3) per-iteration costs 
(4) rate of convergence 
(S) faults oz Jimatations 


(6) advantages over other methcds 


The following sections describe the methods by 


which the parameters « and d are chosen. 


From the Taylor series expansion of f(x) it can be seen 
that not only is -Vf£(x) (if Vf#0) a usable direction but it 
is the locally steepest direction of descent. This first 
order gradient method, though it selects the locally 
steepest direction, possesses no other optimality properties 
such as finding the minimum in the least number of steps. 

(lcelect auitial trial point x =x. x: 

(2)Generate a direction of search d“=-Vé (x*) 

(3) Find the value of &,>0 minimizing f (x*+ xa‘) and obtain 
the next trial point as 

xkti=xheu a‘; 


(4) Set k=k+1 and repeat from step (2) as required. 
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Figure 2.1 Steepest descent search 


Cost Factors: 


(1) Easily programmed and low storage requirement (of order 
n) ; 

(2) Low per-iteration cost; 

(3) Usually possesses a high initial rate of convergence; 

(4) Prone to zig-zagging (see Figure 2.2), decelerating the 
rate of convergence as x* is approached for functions 
which are more than slightly eccentric, thus requiring 


a large number of iterations. 


2-3 Parallel Tangents (PARTAN) Procedure [35] 


PARTAN is a technique of accelerating the steepest 


descents method. 


If f(x) is quadratic and x€E2, then the minimum f (x*) 
can be located in one step, after x! and x? are generated 
from x° by the steepest descent method. This is done by 
searching for the minimum along the ray defined by x2-x®, 
This gives rise to the following algorithm: 
Definition of Algorithm 

(1) Select initial trial point x*=x-1=x9 eX; 
(2)Generate the direction of search d according to the 
following rules: 
d0=~ VE (x°) ; 
at=-Ve(x*) for k=1,3,57Teees3 


qd =x* -xk-3 for K=2 4767) wees 
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Figure 2.2 Zig-zagging of SD along a narrow valley 
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(3) Find the value of 2&,>0 minimizing £(x*+ 0a‘); and obtain 
the next trial point as xktt=x*+u4*; 

(4) Set k=k+1 and: 
(a) for Continued Partan repeat from step(2) as 
required. 
(b) for Iterated Partan continue from step(2) until 
KS 2uit l. Then tre-initialize setting x*=xer-1, k=0 and 


continue from step(2) as required. 


If £(x) is a positive-definite quadratic 
£ (x) =x" Qxta’xtb 

then:[ 22] 

(a) the gradient directions (d°,a1), (d1,d3), 
(d3,d5)-,+s. are orthogonal; 

(b) the paid of) wectors x ,x*7< Re [yey yaipetere 
ane Q-conjugate; 

(c)convergence occurs within 2n-1 iterations. 


(1) low storage requirements (order n); 


(2)low per iteration cost; 


(3) possesses a good initial rate of convergence which tends 


to deteriorate as the minimum iS approached For 


nonquadratics; 


(4)requires a large number of iterations for eccentric 


functions. 
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Figure 2.3 PARTAN search 
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2.4 Fletcher-Reeves Conjugate Gradient Alorithm [14] 


a a ea es we ee et SS ae ee Se Se ee ee eS ee 


This first-order gradient method uses a general form of 
the Gram-Schmidt orthogonalization process to generate a set 
of conjugate descent directions. The method converts 
-V£(x*) into the direction a“ such that it becomes O- 
conjugate with the previously generated d's. These 
directicns are generated by: 

d0=-Vf£ (x9) 

at =-VE (x") +) Bia’ —1 k>1 
ue its are chosen ? Satisfy the condition that a“ should 
be conjugate to all previous directions. The general 
formula for B' can be established [22] as 

=< VE (z") , VE (x \>/< VE (x*=1) ,VE(x1)> i= 


et i#k 


Practical experience has shown that the performance of 
the method is improved if it is restarted from time to time 
due to contamination of B dk-1 by computational error. 
However, the restarting scheme often turns out to be problem 
dependent. Still, restarting every n+1 iterations 2S 
generally the best scheme [14,15,22]. 

Properties: 

If £(x)is a positive definite quadratic 
£ (x) =x" Qxta’x+tb 

then: 


(a)the generated directions are Q-conjugate; 
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(b) convergence occurs within n iterations; 
However, due to round-off in the arithmetic process and 
inaccuracies in determining x, it usually requires more than 
h iterations and the procedure must be restarted according 
to some formula. 
Definition of Algorithm 
(1) Select initial trial point x*=xo <x; 
(2)Generate a direction of search according to: 
d0=-Vf (x0) 
d= =VE (en) e608 > Kat, 2, 370 ae 
Where B =<Vé (x") ,V£ (x*) >/<VE (xk-1) , VE (x*-1) >; 

(3) Find the value of «, >0 minimizing £ (x + 0d") and obtain 
the next trial point as 
xkti=xkeq ak; 

(4)Set k=k+1 and repeat from step (2) as required. [If 
k=n+1 then reinitialize setting k=0, and x9=xK+t+1 and 
continue from step(2). 

Cost_Factors 

(1) low storage requirements (order n) 3 

(2) low per iteration costs; 

(3) requires very accurate determination of «; 


(4)sensitive to re-start policy, which is usually probiem 


dependent and for which no set guideline exists. 


This second-order method makes use of the Hessian, 


A(z)? *of vE(x) invcomputing thesdiréctions of -séarchs. A (x) 
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must be positive-definite. 


Consider a second order expansion of f(x) about x: 
£ (x) =f (x*) + (x-xky” Ve (xk) eb xexky a(x) (x- xk) 
If this is differentiated with respect to x, one obtains 
V£é (x) =V£(x*) +A (x*) (x-x*) 
If this is computed at xkt1, then, trying to maximize £(x)- 


f(xk+1) at the xk 


ateration;, one! picks) «x"t£ | so thar 
f (xk+1)=0,. This implies 
x ke1sy* -a-1 (x*) VE (x*) 


Tf £(x) is a quadratic, then a single step of unit length 


along the direction gives the solution. 


The method is very effective for a large class of 
problems. For difficuities and improvements of the method 
see [12]. 

Definition_of Algorithm 
(WWiselect anitial trial point x*=x° x: 
(2)Generate a direction of search 
a‘ =-a-1 (x*) Ve (x*); 
(3) Find the value of &>0 minimizing £ (xk +u,d*) and obtain 
the next trial point as 

xkti=axk+au,a*; 

(4) Set k=k+1 and repeat from step (2) as required. 
Cost_Factors: 
(1) converges to a minimum very rapidly for a large class of 


preblems; 


(2)A-1 may not always exist or may be ill-conditioned; 
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(3)in a non-convex problem £(x) may not decrease when x is 
not near the minimum; 

(4) it is often tedious or impossible to calculate the 
analytic forms for A; 

(5) high storage requirements of order (n2/2); 

(6) high per-iteration cost due to matrix inversion required 


£Or A=2 (2). 


2.6 Variable Metric or Quasi-Newton Methods 


These methods make use of first order gradient 
information, assume that the function is at least locally 
quadratic, and usually attempt to estimate A-1. The methods 
consist of calculating an nxn matrix H and forming a search 
direction d by: 

ak =-4* Ve (xy 


(for later convenience define y<=V£E(xk+1) - Ve (x*)). 
2.6.1 Davidon-Fletcher-Powell (DFP) [7,13] 


In this version H9° is any symmetric positive-definite 
Matrix, usually the identity matrix I. The iterative 
formula for calculating H*+1 is: 

qk+1ay*k+4 «[aar)- (Hy) (Hy)’ 
7 a Net 


yey 


20682 Selt-scaling Variable Metric (SSVM)p . [29] 


In this version H° is any symmetric positive-definite 


matrix, usually the identity matrix I. The iterative 
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formula for calculating Hk+1 is: 


k kK k 
uk+a= afy'a | |H-(Hy) (Hy)” | + «|aa™ 
yTH yTH y yTa 


2.6.3 Revised Variable Metric (RVM) [12] 


es eae ie tS 


in this version H®° is any positive semi-definite 
matrix, usually the identity matrix I. The iterative 


formula for Hk+t1 jis: 


Hk +1=y" 4] (xd-Hy) (“d-H ‘| 


ja In y 


K 


2.6.4 Properties of the Variable Metric Methods [{ 12,13,29] 


me a a ae ae Ss ae See 


For a quadratic these methods reach the minimum point 


x* and, in general, H ccnverges to A-1 in n iterations. 


SSVM has the following properties in addition to those 
of DFP: Invariance to scaling of the function or variables, 
apparently better convergence for large n(>6), and less 
sensitivity to round-off errors and to the accuracy of the 


one-dimensional search for «. H does not converge to A-!, 


Fiacco and McCormick [12] state the following for the 
RVM method: Because H is positive semi-definite, one is 
allowed to use the variable metric method on large problems 
with a structure that makes it convenient or possible to 
compute second partial derivatives for some of the variables 
and not for others. This may be done in the following 


fashion: 
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and the change from HK to Hk+1 is as stated in the updating 
formula for RVM method. This advantage allows H to approach 
A-~1 in fewer than n steps. This method is not as sensitive 


to the accuracy of the search for « as DFP. 


For other schemes cf updating the appoximation H 
see: [21,737,356 4< 

(1) Fast second-order convergence near a minimum and finite 
convergence for a quadratic functional; 

(2) storage of order (n2/2); 

(3) faster convergence for general non-gquadratic functions 
when compared with other quadratically convergent 
methods; 

(4) first-order method so there is no need to calculate 
seccnd-order derivatives; 

(5)the methods are basically stable; 

(6)the methods may become unstable if the Hessian is 
Singular near a solution, but this can usually be 
avoided if the methods are restarted; 

(7)DFP and most of the variable metric methods require very 


accurate determination of &. 
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2.7 Optimization Along a Line 


The one common feature that all the above methods have 
is that they require a search for local minima along each in 
a sequence of directions. Since the methods call for 
humerous linear searches it is very important to do them 


Simply and efficiently. 


The one-variable optimization problem in E! is: Given 
x* and a search direction d* find «, such that 


£(x"+ 40°) =min f(x +0") 


in general, the methods search for «, in some bounded 
interval [L,U] using some strategy to approximate «,or a 
subinterval containing ~,. The upper and dower limits, 1. 
and U, can be found by a number of different methods: 
Quadratic extrapolation and/or interval doubling, using an 
increasing power of the Fibonacci ratio, or simple muitiples 


of the previous upper limit. 


The efficiency of the linear search is very dependent 
on the efficiency of the method used to determine U, which 
in turn is very problem-dependent. All of the above methods 
can be modified slightly to increase their efficiency in 


specific problems. 


2.7.1 Fibonacci_and Golden Searches [12,22] 


Of the sequential search techniques, the theoretically 


most efficient algorithm for fanding the subinterval 
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containing the minimum of a unimodal function is called the 
Fibonacci search and its simplified version is the Golden 
Section method. It is optimal in that it gives the smallest 
Patiosot final: tolinizial interval for a fixed numbeu )-6f 


function evaluations. 


The Golden Section method is a numerical approximation 
to the Fibonacci search and is not optimal in the same 
sense. It ais preferred because the number of function 
evaluations need not be specified in advance and the 
Fibonacci numbers need not be computed or stored in the 


program. 


Definition of Algorithm 

(1) Set the lower bound L=6 
and the initial upper bcund U=0. 
Set counter k=0; 

(2) Find an upper bound on uk using the limiting Fibonacci 
ratio 1.618 = (1+J5)/2 [12]. 
Let ut =u*-14(1.618)" ; 

(3)I£ £(x+0“ a) >£(x+U'-1a) then U“ is the upper bound and go 
to step (4); 
else return to step (2); 

(4jCalculate T,=0.382(U-L) +L and £,-£ (x+T, 4d) 

T,=0.618(U-L)+L and f,=f (x+T, d) 

Wo to step (/) ; 

(5)Calculate T,=0.382(U-L)+L and f, =f (xt+T, 4) 


go to step (7); 
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(6) Calculate 1,=0.618(U-L)+L and f,=f (x+T,d) 
goto tstep® (hy); 

(ARB a fT Serrt=s. 7 PP HT sift 
andsgo 'toustep (6). 
Cpeelecaser UalC et =Lap Lo=is, 


and go to step (3)°. 


dips (T,-T.1 < € then the interval has been reduced to 
the required accuracy and the algorithm is terminated with 


ee ig 7a 


Other termination criteria can also be used. [22]. 


The Golden Search is second only to the Fibonacci 
Search, from which it is derived, as a sequential univariate 
search frocedure that guarantees a prescribed accuracy. 
However, af “é€5 iS cstiall, a large number of ~function 
evaluations are required. Because of this, methods based on 
curve fitting have been found less costly and to give better 
accuracy. The Fibonacci and Golden Searches are recommended 
mainly when the objective function is irregular and 
approximation by a low-order polynomial is likely to be 


inadequate (such as sometimes occurs in SUMT [12]). 
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The minimum along the line can also be approximated by 


curve fitting. This consists of representing the function 
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along a line by a simple poynomial and using its minimum to 
approximate that of the function. This allows us to replace 
a complex function by a simple function whose minimum is 
known analytically. A very precise fit is not needed, 
therefore the function can be approximated by a low-order 
(second or third) polyncmial. This assumes. that the 
function is smooth and unimodal in the interval being 
considered. In most practical applications the selection of 
a smooth enough interval assures us of the validity of this 


assumption. 


Where function gradients are available Davidon's method 
of cubic interpolaticn is the method chosen to locate the 
minimum along a descent direction. This method requires the 
value of the function and its gradient at two points which 
bracket the minimun. 

(1)iR1=<Vird> 
(2) Quadratic extrapolation is used to find an initial 
estimate of the upper limit of &. 
W==2 (£°-£) /R1 

where £7 is’the “value” “of “£ ‘from “the )iprevicus: WFP 

iteration. 

Ancther approximation to the upper limit of “is 

RO= 1/7 (<d,d>) 
Then let “=min(RO, 4 , 16) 


x'=0, ft=f 
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(3) y =x+ad 
(4) T£ £(y)>£(x) or <V£(y) ,d>=R2>0 


Then go to step (6) 


f'=f (y) 
go to step (2) 
(5) h=X%=- of 
Z=3.(£ "=f (yy yk FRAtRZ 
w= (max (z2-R1*R2,0))”% 
d=h (1-(R2+w-z)/ (R2-R1+2W) 
X= f+ ot 

(7) x=x+ 0" 

Calcudate £(x) and Vf£(x) 

(8)If £(x) >min(f',f(y)) then the interpolation 
in order to obtain better accuracy. The 
chesen depends on the sign of <Vf (x) ,d> 
If negative the interval (x",4) is used, 
the interval («,2") is used. 

m 


The value of oo thus obtained is accepted 


approximation to the minimum. 


Repeated interpolation using the cubic 
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is repeated 


subinterval 


if positive 


as a good 


method is 


generally more efficient than the Golden method. The 


efficiency and accuracy is particularly good if the 
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objective efunction: iss “smooth enough to be closely 
approximated by a cubic polynomial. The method can fail if 
the objective function is not smooth enough. If great 
accuracy is required, the computational cost of repeated 
approximation and calculation of gradients will decrease the 


effectiveness of this method. 


If the gradients are very expensive to calculate or 
unavailable then the Golden method of interpolation could be 
used. Davies! algorithn, or Powell's method could 
alternatively be used as the method of minimization 


{see [1,22] for explanation of the latter two methods). 


2.0 Case Studies 
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In this section the results of applying the 
unconstrained minimization techniques presented earlier in 
sections 2.2 to 2.7 to a set of unconstrained problems are 
presented. Most of the test problems are standard in the 
literature and can be found in references 
[12,13,21,22,30,35]. The computer code used in the solution 
of these problems and a guide to its use may be found in 


[49}. 
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In each case the function f(x) is to be minimized, x0 
isatbherstartingypoint and x* /iSe,the) optimum: _pointegwhick 


minimizes £(x) over the feasible set X. 
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Starting Conditions: Each test problem has its own starting 
point x®, For the variable metric algorithm the initial 
metric H®° is taken to be the identity matrix I. 
Stopping Conditions: An algorithm is stopped when all the 
following inequalities are satisfied: 
a d<e 
a2dd< € 
ele a ie Pa peer nl 
Ewas taken to be 10-6 
One Dimensional Search: In the cubic interpolation algorithn 
it is scmetimes desireable to check if the interpolation is 
accurate enough. This can be done by checking if the vector 
product 
<a’ , VE(x'+1)> < «1 
If this inequality is not satisfied the interpolation is 
done again over a refined interval. However, for the 
problems under consideration this waS found to be 
unnecessary and was omitted. The extra calculations and 
function evaluations required resulted in only a small 


decrease in the number of iterations and in a large increase 


in the total time. 


The Golden search was stopped when 
tT 7T.N < € 
€ was taken to be 107-2 for both the unconstrained and 
constrained problems. A smaller value of epsilon resulted 


in fewer iterations but more function evaluations and 
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generally an increase in the total time. 


The problems used to test the programs implementing the 
algorithms described above are given below. Each problem is 
given a number and is referred to by that number. 
TPs) (2-varieble quadratic) 
f(x) ze + Saxe 
x9=(2, 3) 
x*=(0,0) 
£ (x*)=0.0 

TP2: (3-variable quadratic) 
i(Xp=Xet ex et oxe 
xe= (173, 1) 
x*=(0,0,0) 
£ (x*) =0.0 

TP3: (2-variable hyperellifse) 
f(x) =x*+x¢-3%, -2x, +2 
x9=(0, 0) 
x*= (0.90856, 1.0) 
-£ (x*)=-1.04426 

TP4: (Fiacco and McCormick [{11]) 
f(2) (Xe 23) ates 
x0=(10-3,102) 
x*= (3%,0.0)=(1.316074,0.0) 
£ (x*) =0.0 


The function has a stationary point f (x)=9.0 at x=(0,0) 
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2P5: (Beale! s (function. (28 }) 

Ee) (toes (12>) es (2. 255k) (1-x5) 

$42 7625-%, (1=x,57))= 

x0=(1,0.8) 

x*= (3,0.5) 

£ (x*) =0.0 
the’ graph of ‘this ‘function exhibits a narrow valley, with 
its bottom a curve tending towards the line x,=1. 
mp6: (Reosenbrock's Function [12 ]) 

£ (x) =100 (x2=x> i]t (1=x, ig 

xo= (—d22, 1) 

x*=(1,1) 

£(x*) =0.0 
The graph of £(x) is distinguished by a banana-shaped 
Valley, @2LSseborion Curve “given byea, ==x7. 
TP is (eteccher Powell. {12 )) 

F(R TOO (Cx 10a othr )c) exo 

where 

X, =ICOS (ZT) 

X,=rsin (21) 

r= (x2+x2)” 

x9=(-1,0,0) 

x*=(1,0,0) 

£ (x*)=0.0 
This function defines a helical valley, with axis in the < 
directicn. 


TP8: (Powell's Function [12]) 
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£(x) = (x, +10x,) 2+5 (x ,-x,) 2+ (x,-2x,) #4 (10x, -x,)* 

xo=(3,;-1,0;71) 

x*=(0,6,0;-0) 

f£ (x*)=0.0 
This function is -distingtished by a’ singuiarity: of) £5 (=*) 3 
it is net positive definite. It was designed especially for 
testing algorithms with a quadratically converging search. 
TP9s (Wood's Function r11y 

Pave WeO (eeox net (1-5) ) <+90 (nex) = 

Fitae ete (C- x Vet (ia ee) 


#198 (1>x5) (1-x,) 


Rees 1, 17 1) 
£ (x*)=0.0 
This problem has a stationary value of £(x)=7/7.876 at 


X= (-0.9679, 0.9471,-0.9695, 0.9512) 


2.8.3 Experimental Results 
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The results of applying the minimization algorithms to 
the above test problems are given in the tables below. 
These results allow conclusions to be drawn concerning the 


relative effectiveness of the algorithms. 


Each algorithm is applied to each test problem and the 
following three parameters are recorded: 


(1) IT- the number of iterations it took for the method to 


converge. 
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(2) #F- the total number of function evaluations taken by 
the algorithn. 

(3) Time- the time in seconds that it took for an IBM 360/67 
to solve the problen. 

The names of these algorithms are abbreviated as follows: 
SD- Steepest Descent 

GNR- Generalized Newton-Raphson 

TPTAN- Iterated PARTAN 

CPTAN- Continued PARTAN 

CG- Conjugate Gradients 

DFP- Davidon-Fletcher-Powell algorithm 

SSVM- Self-Scaling Variable Metric algorithm 


RVM- Revised Variable Metric 
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peee (eer, 326 | He § 6t2 (te Th ae 
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Table 2.1 Cubic interpolation 
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* =- METHOD DID NOT CONVERGE TO CORRECT SOLUTION 


Table 2.2 Golden Search 
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From the tables: it can be seen that the restiits ere 
basically problem-dependent, but that certain general 


conclusions can be drawn. 


The three variable metric algorithms and the conjugate 
gradient algorithm solved the two quadratic test problems in 
n iterations. (The table of results shows n+1 because that 
is the minimum allowed by the program though the methods 
actually converged earlier.) Both of the forms of PARTAN 
required 2n feerat ons: for these problems. These results 
agree with the theoretical expectations for these algorithms 
for quadratic functions. Steepest descent required many 
more iterations but because of its simpler logic which 
requires fewer calculations, the total time required for 


solution is only 40% higher. 


It is in the sclution of the non-quadratic problems 
that the greater eficiency of the accelerator techniques 
nce apparent. From Tables 2.1 and 2.2 2% can be seen 
that in general SD is less efficient in terms of number of 
iteraticns and total time by one order of magnitude. It is 


also apparent that PARTAN is much less efficient than the 


other techniques. 
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algorithms are approximately equal in performance in terms 
cf number of iterations and time. C-G generally takes more 
iteraticns but less time because of its simpler algorithn. 
The variable metrics are generally equal in performance but 


each has at least one problem in which it excels. 


For test problems 4 tc 9 the performance of all the 
algorithms was found to he very sensitive to the mode of 
pailenentasida of the one-dimensional search algorithms. 
Their performance can be varied by factors of 2 or 3 both 
adversely and favourably. Ali numbers in the tables were 
obtained from the same program, that implementation being 


the one that gave the best resuits in general. 
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Two methods cf one-dimensional search were used: Golden 


Search and Davidon's cubic interpolaticn. 


From Table 2.2 it can be seen that for the same number 
of iterations the Golden method takes 25% to 100% longer 
with 4 to 6 times as many function evaluations as the cubic 
method. The only saving here is that only the final 
function evaluation of each iteration requires the 
calculation of the gradient. For several of the examples an 
e of 10-2 is not accurate enough to achieve satisfactory 
results: but decreasing ¢€ increases the number of function 
evaluations per iteration so that while there would be a 


decrease in the number of iterations there would not be a 
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correspcending decrease in time. 


in general, one can say that the cubic interpolation 
algorithm gives more rapid and dependable convergence than 
the Golden Search. In addition, the performance of the 
Golden Search is very dependent upon the individual problen 
and the choice of ¢ for which no a priori guidelines exist. 
The main advantage of the Golden Search is that it does not 
require the calculation of gradients and its accuracy is not 
dependent upon the sometime questionable accuracy of a_ low- 


order pcelynomial fit. 
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Convergence of the Algorithms 


Since the SD method has a high initial rate of 
convergence, and since the convergence of the variable 
metric techniques is sometimes hampered by contamination of 
the metric due to poor initial search directions, it was 
decided that trying a few initial iterations of SD at the 
beginning of each search might improve the rate of 
neon tee Tables 2.3 and 2.4 present the results of 


trying three initial iteraticns of SD. 


It was found through experiment that three iterations 
was the best number. Less made little difference and more 
often hampered convergence though sometimes convergence was 


improved. 
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From the tables one can see that SSVM and conjugate 
gradients are relatively unaffected; both methods of PARTAN 
and RVM are both adversely and beneficially affected. 
However, DFP generally benefits from the initial iterations 


of SD searches. 


In conclusion, from the evidence one gains the 
tupeession “that “thes veffiect! of “initial aAterations of SD 
searches is problem dependent but that the DFP algorithm 


usually benefits from it. 
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3. Constrained Optimization Strategies 

In this chapter some of the methods of solving 
constrained minimizaticn problems are reviewed. Emphasis is 
placed cn the SUMT algorithm of Fidcco and McCormick and on 


Rosen's Gradient Projection Algorithn. 


In the real world one cannot always solve a 
minimization problem without running into constraints. In 
dynamic-system design, for example, there are usually limits 
on loads, velocities, accelerations, etc., which result in 
constraints being placed on the solutions. Unconstrained 
methods cannot cope directly with these limitations and so 
different techniques Or reformulations of existing 


techniques must be used to solve these problems. 
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The MP problem involves the minimization of an 
objective function, this function being subject to 
inequalities or equalities called constraints. In general, 
the objective function and constraints may be nonlinear 
functions in all or some of the variables. The special case 


of no ccnstraints was handled earlier in Chapter 2. 


Any set of variables which satisfies all the 


constraints must be a feasible solution, and the feasible 


ay ae LT ‘i = 
if 


tm f i 
gale fer ae ition gil : he ii Power) hie a 
(i) olan ty eb yo estate be oO 
je bee to deta tiga | pay Lonmin ys gil tex0.28 win ae se nes 


nay = 
7 2 


invshaniyll AO (* oot ont Mite aa ie Mat 
ere =Viiw iy meu aee ae i@ou jaa?” on ne — 


‘7 
4 


a hep teese sy wh gh paifands. tur i928 wetloxg wilsatte ft ke. 
5 L FRLcan wae eral oy (GRewa: to8 sin PeeL “ae ‘a 
i URGeT P98 LN ate iiektarsiowss «3ehes Lpolew vito k ~ 10 
hee tas saquié sea EnioS) oft de Cetaly ofis7 ota bastene > 

ae Stelonrhack Snel Biiv Qiese Ste ages sents ohue-en 


ia fb2aBe T9 AALS? BarwIi4 i sav cahoe ?t *19 2932 7 


netiadga wamis Wel @ vi tea at Suva weep ieavee fi 
Wit 


‘— "1% ) ie hs em@-oaspesaues., bogétogndsemeds Fi 


cc 
7m 2 GOloeetoaase ih: mrhivalt aa«ftiey " aa? 
vg +[ePatie Garr pork? ging elds ,eOltoaw.§ ave sos 
‘ ; a 


~levoney af Atebes tebe. bndfeo: we) *eheepe Goi cuels? 


a r 


at. 
shoal 
on) 
Tom 


te ee 2. ae tase | Ona Spttouds pai. 
| Pees 
: 


easy lascejya ae meiveccey OF 1 eens ie tte a atoltoie 
Ca gehynat - re cham) Pade deed atnkesoan ‘+0 
adh ba pores aati Wiodaly we i, 


ota, ee tbe aia ee rie waa 
~~ a ne | | 


39 


solution which optimizes the objective function is the 
optimal soluticn to the ME problen. 
The general form of the MP problem is: 
"Minimize f (x) 
Subject to g/-(x)s0 1=Vyreee ell 
h; (x) =6 Tt 4 eee 
The scalar objective function f(x) is a function of the n- 


dimensicnal vector x. 


There are several special classes of the MP problem. 
aiie(x), o(x), and h(x) are all linear then “one hasS the 
linear programming (LP) problem. If any of £(x), g(x) or 
h(x) are nonlinear then one has the non-linear programming 
(NLP) problem. The quadratic programming (QP) problem, a 
special case of the NLP problem, has a guadratic objective 
function with linear constraints. There are other special 
types such as the geometric programming problem, integer 
programming problem and stochastic programming problem: They 
are not considered here. There exist special technigues for 
the solution of the QP and LP problems but they are not 


involved in this study. 
3.2 The Kuhn-Tucker_ Conditions [20,44] 


Consider the NLF problem in the form stated above. at 
x* is an optimal soluticn to the NLP problem, then the Kuhn- 
Tucker (K-T) conditions (under reasonable assumptions) are 


necessarily satisfied. Generally the K-T conditions state 
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that given the point x*, if one moves in any direction da 
distance xj) dif (&>0);, as long as one remains in ‘the 
feasible region the objective function will not decrease. 
In addition "The K-T theorem is intimately related to the 
Lagrange multiplier theorem, which it generalizes,... either 


theorem can be used ta obtain the other."f{ 20] 


The K-T conditions for the above NLP problem are: 
If x* is an optimal solution and a suitable "constraint 
qualification" holds at x* then the following condiions also 
hold: 
(1) x* is Peeet pies 
(2)there exist multipliers 4;20, L=1 piers eg hl 
and unconstrained Ce ee Ave L=Mt 1 po<erG 
euch that: £ (x4) + ) ign (x4) =0; 
and “ 
(3) ig; (x*) =0 ret hrc er ral 
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The K-T conditions are necessary but not generally 
sufficient criteria) for ddentityingcan-optanaly pointig sas 
such they are important in theoretical developments and 
serve as a basis for many computational algorithms. They 
also have the added attractiom that they constitute a first 
order necessary condition for a minimum which can be tested 


numerically. 


3.3 The Constrained Problem Transformed to an Unconstrained 


ae ee re a ae ee a a ee a ee ee eee OOO OOO OO CO OO eee eee ee 


seGo Setagy Ww aioe ad? Go¥ Wooten aires 
ees To i bis s) Ae fab antes Leatsqo a aa 
i 45) ion ps Sap Ree one wuts ce 3p dhtow 


T puee ge el eo | uy @1¢ pay sive Foixe one 
J 2 nie 
ead s et OOP k yh eam ‘fine TAd) wr ranodsew Oho | 
E2= Me) gach © ere Ay :t04? t9084 Se 
—"? i _ 
hes 
- 
er ™ de (*a) ote A (E)/= - 
l - 
OSS ateuthi tag "Hii rante io" 4 ApeeMmsan T-3 wae a ad 
A Peet = Ae nee we Oe Me 5, stein ae 7 
fr’ £2 HERTS LARD 1) Ae ys eae >. £267 guns - 63h vads doe’ 


Geil Sky Sek calityL esenpaiaes val an | PA ike ae s .*6 oe 


3s @ & 67712 %*t in! 


we F 


Th 


The solution to the general NLP problem can often be 
found by replacing it with a sequence of unconstrained 
problems, each containing the constraints implicitly, whose 
limiting solution is the NLP problem solution. This 
unconstrained optimization (indirect) approach is motivated 
by several practical considerations. These are: 

(1) The unconstrained formulations of the constrained 
preblems are often simpler to apply or use; 

(2)Many more efficient algorithms have been developed to 
handle unconstrained problems; 

(3) The unconstrained algorithms have been more extensively 
studied and further develcped; 

(4)One avoids the necessity of having to cope individually 


with the boundaries of the feasible region. 


The disadvantage of these methods is that they are not 
as efficient fer some problems as the direct methods since 
the structure of the original problem is obscured. This is 
especially true when one has linear constraints or a linear 
objective function. The resulting solutions may also only 
be poor approximations to the correct solution instead of 


exact. 


One way to modify the objective function in order to 
incorporate the influence cf constraints is to add something 


to it which would bend the search direction into the 
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feasible region. In penalty (exterior point) methods one 
attempts to create an infinite penalty for departing far 
from the feasible region, and in barrier (interior point) 
methods one attempts to set up a barrier against leaving the 


feasible region. 


The following are forms of the Penalty-Barrier methods: 


man ¥ (X,o)=£ (x) +27 TP (x) 
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koe 
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UM VE (x, .) =P (x) eB (x) 
lin r,=0+ 


K+ 0o 


With Tespect to.all x e {x}" 
3.4.1 Penalty Constrained (Exterior Point) Methods [12,20] 


By construction, the penalty function P(x) is ideally 
zero for x €X, positive for x 4X, and wsually forces £ (2,2) 
toward +co for r~0O and mers The typical penalty function 
strategy is to solve a sequence of penalty-constrained 


problems with r given successively smaller positive values. 


fv pical ly the penalty functions are defined for 
inequality constraints as; 
P, (x)= 0 for 9g, (x)ce 
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"Min £(x) Such that a<x<b" 
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Figure 3.1 Penalty method 
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Figure 3.2 Barrier Method 
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This can be formulated as 
Pi (x) =( (9; (x) + ilo; (x) I) 72) 2 
or as 
P) (x) =Max2 (0,9; (x)} 
For equality constraints 
es (x) =h2 (x) 
Thus, given the NLP problen 
"Min f(x) 
subject to g (x) <0 
h (x) =0" 
mM Y 
F(x,)=£(x) +r () gi (x) +> bi (O) 
Ge L=m+i 
The advantages of the penalty-constrained method are 
that the initial point x°® does not have to lie within the 
feasible region and that it handles both equality and 


inequality constraints equally well. 


The disadvantages of the method are: 

(1)it does not produce intermediate solutions that are 
feasible; 

(2)the problem becomes ill-conditioned as the penalty 
coefficient increases; 

(3)it is often difficult to minimize as there is a lack in 
the orders of differentiability in x of the function at 
any boundary point cf the feasible region; 

(4) numerical errors in the penalty terms become Significant 


for large penalty coefficients. 
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3.4.2 Barrier Constrained (Interior Point) Methods [12,20] 


For xeéint{xX} the barrier function B(x) is continuous, 
tends tcwards -o as any component of g(x) tends to 0O-, and 
ee ye (=) ranges on (-~,0 J. AS with the penalty 
constrained problem, the strategy calls for solving a 
sequence of problems, in which the parameter r is assigned 


successively smaller positive values, 
Typically B(x) is defined as 


(1) B(x) = 2,972 (x) 


Lat 


or (2) B(x)=)_Un(-g; (x) 


v= | 

Thus, the NLP problem 
"Min £(x) 
subject to g(x) 2o" 
becomes 


F (x,xr) =f (x) -rB (x) 


The advantages of the barrier constrained method are: 
(1)it produces intermediate feasible sclutions; 


(2)it gives a final solution of arbitrary accuracy. 


The disadvantages of the method are: 
(1)an initial point in the feasible region is required; 
(2)the rapid change of B(x,r) in the vicinity of the 
boundary complicates the one-dimensional optimization 


problem which most techniques require; 
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(3)at the boundary B(x,r) gets very large and r,~-0+ so the 
function is very sensitive to variable Changes and 
round-off error is introduced; 


(4)it does not handle equality constraints. 


3.4.3 Combined Penalty-Barrier Problems [12] 


There are many circumstances under which Lt is 
desireable to use a combined penalty-barrier algorithm. If 
one wishes to take advantage of the superior behavior of the 
interior point methods but also has equality constraints, 
then the mixed algorithm must be used. It is also necessary 
to use it in penalty constrained problems when continual 
Satisfaction of certain constraints is necessary (such as 
when one of the constraints or the objective function 


eontainsva logarithnic*tern).< 


This problem takes the form 
WMan F(x,©C)=£ (x)+r- tP (x) -5rsB (x) " 
In this method the initial point must lie within the 


feasible region defined by the constaints contained in B(x). 


3.5 The_ Sequential Unconstrained Minimization Technigue 


SS Se — 


{sumMT) [12] 


SUMT was developed by Fiacco and McCormick [9,10,11,12] 
based on an earlier technique proposed by C.W. Carrol. Lt 
is a mixed penalty-barrier technique which deals efficiently 


with non-linear objective functions and both equality and 
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inequality constraints. It generates a sequence of 


solutions which converges to an extremum. 


The convexity or concavity of the functions involved is 
not critical to the performance of the algorithm; however if 
certain conditions are not satisfied it may converge to a 


local instead of a global extremun. 


The following conditions are necessary for a local 

solution: 

(1)the interior feasible region is not empty and an initial 
interior feasible point x° can be found; 

(2)the functions £(x) and g(x) are twice continuously 
differentiable; 

(3)£ (x) is bounded below on the feasible region; 

(4)all local minima are located at finite points; 

(5)if ain addition £(x) and g(x) are convex and at least one 
is strictly convex, then the method wiil converge to 


the global. ninimun. 


Given the problem 
"Min f(x) 
subject to g (x) <0 
h (x) =O" 
then the objective function is defined as 


m ¥ 
F(x,r) =£(x)-r ) In fg (x) ) #r-1 ) -h2 (x) 
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(1)Given xe Xeand 150 set k=1; 
(2)Minimize (using an unconstrained minimization technique) : 
m Y 
F(x ,r)=f (x) -r, ) Ln (g; (x)) #r¢4) (he (x) ; 
v= Lz M+( 
(3) Set r,.., =r,/fce where c>1 is scme constant and k=k+1 
(4) Test for convergence: 
oP tae beeen Cer 
WE(X,Te4, ) Ml 
go to step 23; else stop. 

Comments 

(1)Other criteria, such as testing for the difference in 
magnitude between successive values of x, can be used 
for termination; 

(2)The selection of r and c are often critical to the rate 
of convergence. Generally r should not be chosen too 
small and c should not be chosen to be too large or too 
small as otherwise convergence difficulties may be 


encountered. For more information on the choice of fr 


and c see {22}. 
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In this section a direct method for the solving of the 
constrained minimization problem is discussed. This method 


generates the search direction by projecting the negative 
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gradient - VE (x) onto a subset of the constraints binding at 


Xe 


Given gq linearly independent n-vectors €, ez 4e00728G iD 


49 © . 
E , these n-vectors span a gq-dimensional subspace denoted by 


Q. If g=n, then Q coincides with E”. 


The equation e'x=c Gefines a hyperplane H. in E or 


manifold of -dimension n-1. The intersection H, to Hg or q 
of these hyperplanes creates an n-q dimensional subspace, 
denoted by Q. Thus Q is the set of all x with 


i 
(2, ¢é €,) x= Ei x=0 
-| ¢ gqeeeeg a a , q 


Each of the vectors e, to eq is perpendicular to Q, 
therefore Q and Q are mutuaily orthogonal. Together Q and Q 


span the space E”. 


Define the gxn matrix P [33] as constructed by the rule 


A nO Ee 
B=E 4 (248g)? 


The matrix Py is a projection matrix which takes any vector 
in E” into Q [33]. The nxn matrix P, is now defined by 

Pp, =I-E.(E,E,)—1E 

En Chay ae 


: : . . = A 
Thus Py is a projection matrix which takes any vector in E 


into the intersection Q [33]. 


Tt can also be shown that 
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in Figure 3.3 the projection z is thus obtained as 


Z=Pq Xe 


Figure 3.3 The projection z=Pqx 
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Let the objective function f(x) be convex with 
continuous second partial derivatives in the feasible region 
X. The vector x is constrained by m linear inequalities of 


the forn 


Where) A=(a) 47-66 aac] 
The constraints are normalized so that 
(aa ea OS Pa 
These constraints restrict the solution to m closed half- 
spaces. Their intersection is, in general, a convex 
polytope and is called the feasible region X. 
Let A,=[4, an. 22 2Gq | Where a-=(a)7 (a2 n0<+ Gnu) 
bg = (bi pDayee+ bq) 
and Q=[H,/ H, «+. Hy] 


where Q is the set of hyperplanes active at x. 


if the point x lies: at the intersection of q tinearly 
independent hyperplanes H,,H,,---,Hg which are defined by 
the q linearly independent vectorS 4, ,dzreee Aq and the 
components Di rbreeeedy (q<m), then the gq equations 
A, x-b,=0 
determine the points that lie on the intersection of 


H, giljee sey H Given Aor the nxn symmetric projection matrix 


a 


TT {i 
Py is given by Pa =I-A,g (Ag Ag) ~*Ag 


vu 


Thus P, will take any vector in E into Q. 
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Figure 3.4 Gradient projection 
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Normalize the linear ccnstraints 
fest sthe feasibility» ‘of the «initials port s#x%and 
determine the constraints active there 
If some of the constraints are active then select those 
that form a linearly independent set and go to (5); 
Else set Fy equal to the identity matrix and go to (7) 
Calculate (A; 0,122, 
Calculate the projection matrix P, 
Calculate the directicn of search d (see 3.6.4) 
Calculate the vector R by 

R=T, A (-VE) 
and let R, =the maximum element of R then 
Tf fladllse€ or IP, (-VE) <€ 
Go to (10); 
Else go to (11) 
If each element of R is <0 then a minimum has been 
found. If one or more elements of R are >0, then the 
objective function can be further minimized by dropping 
the hyperplane corresponding to R, and returning to (5) 
{33}. If none of the above go to (11) 
Let Beta equal the maximum of the sums of the absolute 
values of the elements of the rows of T, 
If R >Beta drop H, from Q 
Normalize d@ and determine the maximum stepsize, tu, that 
can be taken along d without causing any of the 


constraints to be violated 
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(13) Let y=xt od 
LELAdy VE(yoa20 sets x=yrand govtor (14); 
Else interpolate to find the y that minimizes f along d 
and «set x=y, then go to (5) 

(14) Add to Q the hyperplane which corresponds to the 


stepsize %, and go to (5) 


3.6.3 Calculation Requirements_of Gradient Projection [33] 
(1) Calculation of (Bg Ag )—? by the use of recursion 
relations: 
(a) Dropping a hyperplane Hg: 
Assume that Grea! is known and Ag =[ a, a.ee-. agi J 


Suppose a Square matrix C is partitioned so that 
g p 


G= \C, Ce 
C., Cy 
where C,,C, are square matrices and C,/C, are 


rectangular matrices. Then 
C-1=B= |B, aD 
Bene ch 
It can be shown that [33]: 
G2 P Be nD 
=) 
Applying this to determining (Ag ae , let 
Ls 
B-Ga (ta) 
and therefore B, 2B,,Ba, and B, are known. In 
Partacular, B, 1s a (q-1)x(q-1) symmetric matrix, B, a 


_ 
(q-1) column vector, B,=B,, and B, a scalar. 


ff a hyperplane 4, (l<q) is to be dropped, the 
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relation still applies if the i-th and q-th rows and 


columns of B are interchanged before it is applied. 


(b) Adding a hyperplane Hg = 
Let - i a ec toe 
and D-1=B= |B, ites 
Ee! 38, 
[hen p= Dest 05 tp p-420, 0, 
BPSD S De 
Ba==D7nU 2072 
ye a 
If D=C, Ca 
; i 
and D, =C 2 coe 
M4 1 
then D,=Cg-1 ag =D3 
aie 
and D, =a, a =1 
giving D, =aq Pq-t a g= IP a a gill®=laq- q-! F (2 


This gives: 
B, =D-14+Dr1FF° 
B, =-D51F=53,” 
B, =D; 
where Epa De 
Thus allowing one to construct 
B= (Ag Ag )—2 
(Cle) Shep sacs teats (Ag Ag) ~2 is calculated it is done by 
Gaussian Elimination. 
(2)Calculation of the projection matrix Py? 
(a)The first time Pq is calculated and everytime a 


hyperplane is dropped 
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Be=i=A, (hah, )- 1k, 
(b) The recursion relaticn for the addition of a hyperplane 
ts. [23/3 
Part =Pq ~Fq Age, ages Py Vase, Pg Ag+! 
(c) Dropping a hyperplane from Q: 
e)5 The point x* is a constrained global minimum 
GE 2 (z)etLe 
Sea e (xe) —0 (It Pe VE (x*) I <6) 
and (Ag AG) 723, Vf (x*) =R<0 
If one or more components of R are positive, this 
indicates that the objective function can be 
decreased by dropping the hyperplane corresponding 


the the largest positive component of R, R,, 
Da ae (WP, VE (x) > &, it may still be desirable to 
drop a hyperplane. Comparing the lengths of 
see ve(y) coals ee Vfi(y) in Figure 3.5 one can _ see 
that more is tc be gained by moving along 
2, , 


be dropped from Q. To determine if this should be 


Vf£(y) than along -P, Vf(y), therefore H, should 


done one calculates the sums of the absolute 
values of the elements of the rows of (Eek)! and 
sets Rk equal to the maximum of these sums. If 


Mm 


Ry>?, then drop H,, from Q. 
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Dropping a hyperplane 
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(3) Calculation of the maximum stepsize t 
Given y=x+td 
we wish to find the maximum value of t% for which all 
Constraints are’ satisfied. At the point where this 
line intersects the hyperplane H. 


y=x+ tid 


substituting for y and solving for 7; gives 
t= (b;-alx)/aid 
(Tis -calcuhated for all hyperplanes not in -Q; The 
minimum positive value of t,is the maximum stepsize 
which can be taken. 
(4) Interpolaticn Routines: 
If the stepsize is infinite then there is no constraint 
in the direction of d and cubic interpolation is used 
to find the minimum cf f along d. 
ie the stepsize is finite then repeated linear 
interpolation is used until 
WaVE(y)I<é 
where y=x+txd («found by interpolation) 


(cubic interpolation could be used here instead) 
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(1) Pure gradient projecticn: 
This is the method originally proposed by Rosen and 
consists of projecting the negative gradient onto the 


active constraint set. The method possesses the poor 
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convergence characteristics of the Steepest Descents 
algcrithm and has only simplicity to recommend it. The 
direction of search is given by 
d=-P, Vf (x) 
(2) PARTAN Acceleration: 
Here the PARTAN acceleration method is adopted in order 
to improve the convergence of Rosen's algorithn. Only 
the Continued PARTAN algorithm is used to determine d. 
In it, every time a hyperplane is either added or 
drcepped the acceleration routine is re-initialized. 
See the discussion of PARTAN in 2.3. 
(3) Conjugate-Gradient Projection Acceleration: 
Here the Fletcher-Reeves C-G algorithm is extended to 
Rosen's method. This extension was suggested by 
Goldfarb [16] and was implemented practically in this 
study. The directions of search are determined by 
a! =-p, Ve (x') +P, V£(x') 12a" Vp, VE (x! IL 2 
while the initial search direction is given by 
doo ae, Vf (x°) 
Every time a hyperplane is either added of dropped the 
acceleration routine must be re-initialized. eye aS 
also re-initialized after n iterations. As long as the 
constraint matrix remains unchanged, this method 
possesses the excellent convergence characteristics of 
the unconstrained C-G algorithm. 
(4) Davidon-Fletcher-Powell Variable Metrics: 


Here the DFP algorithm as developed by Goldfarb [16] is 
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extended to Rosen's method. 
The directions of search are given by 

=-HVE (x) 
in the first iteration H® is set equal to Pg « 
Whenever a hyperplane is dropped H is updated by 

7 ae 

He, Ha #Ps-) Aq Gq P52) /aq Poan 
and whenever a hyperplane is added H is updated by 

Has: Hq Hq @quitqes Hq /Aq.(Ha agus 
When the constraint set remains constant H is updated 


by 


where z=td 
Y= Vte je x”) 
B=zz' /z'y 


c=", yy" H, /y H;Y 


This extension has been found to improve greatly 


the convergence characteristics of Rosen's algorithn. 


In this section the results of applying the metheds of 
constrained minimization discussed in Chapter 3 to a set of 
constrained problems are presented. No claim of 
completeness or special difficulty is made for these 
problems; they are presented, in the main, as typical 


examples. 
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in each case the functional £(x) is to be minimized, x9 
a8 the starting point and x* is the optimum point” Which 
minimizes f(x) over the feasible set X. 
Starting Conditions: Each test problem has its own starting 
pornt —'x?. For the methceds using the variable metric 
algorithms, the initial matrix H° is taken to be the 
TOeHTICy Tatrix 1. 
Stopping Conditions: These vary with each method. See the 
discussion of experimental results for each problem. The 


computer code used in the solution of these problems and a 


guide tc its use may be found in {49,50}. 


The problems used to test the programs implementing the 
methods in Chapter 3 are given below. Each problem is given 


a number and is referred to by that number. 
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TX PX¢ 
Gq (X) =" Xe 
(aepas Oe 1e7e, wpa) (God A160 7065 705707 70 201073.0 10) 
¥e='(10207-6.57.0. 174202120) 
x¥=(4.0,1.0,0.5,0.5,3.001) 
£ (x*) =0.540 
This NIP problem gives the solution to a third-order 
nonlinear system whose transfer function is 
G(s,P)=1/[ (sta) (s2#+2dwstw2) ] 
p=(a,d,w) >0 
The requirement is to maximize the gain K=1/x, {or minimize 


xX,) and still retain the stability of the systen. 
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x* =(0.947,1.096,0.2245,-0.0864) 

£(2*)==1,058 
This problem has convex constraints but a non-convex 
objective function. Examination shows that there are, in 
addition to the origin, two other points at which the Kuhn- 
Tucker conditions are satisfied: 
X=1(0 09 - 2/37 9723) 
£ (x) =-0.222 
and 
x= (0,0,0.707,-0. 754) 


£ (x) =-0.533 


13: Colville's Test Froblem One [22] 


ie 


£(X)=- 15x, ~27x,-36x, - 18x, - 12x, 
qoOxe 40x x) - 20%, x2 tO4x) x, - 20K xe 
+39x2-12x (x, ~62x, x, +64x, x, 
+10x2- 12x, xX, 20K, x +39x2-40x xX +30xK2 
TAKS +E XS t1OxS + OXI + 2x3 


subject to Ax<b where 


16 -2 0 -1 0 
0 2 0 -0.4 -2 
3.5 pie b=2Z 0 0 
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b =({40 (270.2574 ,4,1,40,60,-5,-1,0,0,0;,0,0) 
x9=(0,0,0,0,0) 
x¥= (0.3,0.33347,0.4,0.42831,0.22396) 


£ (x*) =-32634868 


3.9 Numerical Studies 


ee se eS a SS ie Ss 


The practical results of solving the constrained 
problem by transforming it into an unconstrained problem are 
daseussed “in sections 3.9.1 to 3.9.3 below. in section 
3.9.4 Rosen's gradient projection algorithm is discussed and 
its performance on TP13 is compared to a solution of the 
Same problem by the SUMT algorithm. The conclusions drawn 


from the numerical studies are presented in section 3.9.5. 
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The penalty method was applied to the solution of TP10, 
TP11, and TP12. DFP, SSVM, RVM, CG and PARTAN were applied 
to the solution of all three problems. The form of the 


method chosen was 
MM an B(x, Cp=f (x) FES ee (x) 


Where P (x)=) -max2 (0,g; (x)}" 
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Table 3.1 Fenalty solutions r=0.2 
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Table 3.2 Penalty solutions r=0.02 
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Table 3.3 Penalty solutions r=0.0902 
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Generally the solutions were inaccurate with the best 
solutions having errors in x of approximately 10%. For 
several cf the problems the solutions were so inaccurate as 
to be completely unacceptable while others converged to 
local minima. As expected, a decrease in the value of fr 
usually resulted in an increase in both accuracy and the 
number of iterations required to solve the problem. The 
variable metric algorithms experienced some difficulty in 
obtaining solutions and in some problems the algorithm 
itself failed resulting in no solution. This, however, 
could be corrected by restarting the searches at suitable 


intervals. 


3.9.2 The Barrier Method 


The barrier method was applied to the solution of TP10, 
TP11, and TP12. DFP, SSVM, RVM, CG and PARTAN were applied 
to the solution of all three problems. The form of the 


method chosen was 
"Min F(x,r) =f (Xx) -rB (x) 


Al 
Where B(x) =) ,1/9; (x)" 
=| 


for several different values of r. The results are given in 
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Table 3.4 Barrier solutions r=10-2 
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Table 3.5 Barrier solutions r=10-@ 
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Comments 


As with the penalty method, the solutions were 
inaccurate. Here too accuracy and number of iterations 
generally increased as r decreased in value. It can be seen 
that FARTAN is completely unsuitable for use as a 
minimization strategy for barrier problems solved in the 
Manner presented here as the function does not sufficiently 
approximate a quadratic. Again the variable metric 
algorithms experienced some difficulty in finding solutions 
but not as much as with the exterior penalty method. Here 
also this difficulty could be corrected by restarting the 


searches at suitable intervals. 


3.9.3 The SUMT Algorithm 


SS a a Swe SS we a Se SS SS 


The SUMT’ algorithm was applied to the solution of TP10, 
TP11, and TP12. DFP, SSVM, RVM, CG and PARTAN were applied 
to the solution of all three problems. The form of the 


method chosen was 
“Min F(x,1r) =f (x) -IrB (x) 


tal 
Where B(x) =) Ln (g; (x))" 


and r°=1 and c=5. (The algorithm implemented here for this 
study was the basic SUMT algorithm. No attempt was made to 
check fcr and include only active constraints as the logic 
required both for the main program and the function gradient 


Subroutine is very complex and the results were not deemed 
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worth the effort.) 


In addition the following were also studied: 

(1)the relative effectiveness of the Golden and cubic 
linear search techniques; 

(2)the advantage gained when using the variable metric 
algcrithms of setting H9 equal to the H generated by 
the previous unccnstrained minimization instead of 
restarting the minimization by setting H® equal to the 


identity matrix. 


The results of these studies are given in Tables 3.7 to 
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Table 3.7 SUMT solutions: Cubic interpolation 
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Table 3.8 SUMT (cubic): Effect of restarting 
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* ~ METHOD DID NOT CONVERGE TO CORRECT SOLUTION 


Table 3.9 SUMT solutions: Golden Search 


iS ry a | ee Ge ee a (a as (mem reer magn [iw ne ener aman GER P ea 
{ RESTARTED | NOT(RESTARTED | 

{TEST | PROB {| DFP [ Sov" | RVM [| DFP | SSVE  — RvVe | 

ee ee ee 

| (eect hGssteeolee hy Teo. | oo [138 (9-121 


| 
{TP10] #F {| 3349 {13909 | 7349 | 1744 | 2562 | 4009 | 
{ (TIMES. 1691(35.205)/47.1131 4.965] 6.735) 93.872| 


f---—+--—- + —----4 -_ -_--- ------------ +--+ + --+4 
fete) 19-6 So es ie Si | 
(TP11] #F | 600 | €55 | 599 { 450 | 528 | 420 | 
| [TIME] 6.498] 0.548] 0.495] 0.437] 0.496] 0.4131 
f-—------------$ ------ + ------ ------+—_—___ +H 
| Pp tyap 17es| “Sit (ae S3e( woody 70 (mete 
(TP12] #F | 4197 | 984 | 3090 {| 5007 { 1345 | 1695 | 
[TIME] 3.591{ 0.966] 2.463] 4.661] 1.357] 1.577] 
Es ee 22) ee eee eee ee ee | 


* - METHOD DID NOT CONVERGE TO CORRECT SOLUTION 


Table 3.10 SUMT(Golden): Effect of restarting 
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Table 3.11 Test problem 13 
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Experiments were also conducted in measuring the 
effects of varying the values of €,«' used to determine the 
convergence of the SUMT algorithm and in varying the value 
of c used. A standard of €,€'=10-* was employed. Choosing 
a smaller value generally resulted in an increase in 
iterations with very little increase in accuracy. A larger 
value sometime resulted in poor solutions so 10-¢* was 
selected as the best value to used on the problems being 
studied. A value of c=5 was chosen as giving the best rate 
of convergence with the fewest difficulties, KR choice™~ of 
c=3 resulted in no improvements and poorer performance for 
some examples. A larger c created difficulties during the 
initial SUMT iterations. When the minimum was sufficiently 
approached a strategy of increasing c to 10 was employed and 


this gave improved results with no detrimental effects. 


As can be seen from Tables 3.7 and 3.9 the Golden 
Search has few advantages to recommend it over the cubic 
search. The conjugate gradients algorithm displayed 
somewhat superior performance when employing the Golden 
Search but for the other algorithms even though there was 
usually a decrease in the total number of iterations, the 


number of function evaluaticns was prohibitively large. 


The cubic interpolation method as implemented presented 


no difficulties during the linear search. Some difficulty 
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is often experienced with the cubic interpolation algorithm 
when applied to barrier constrained problems. This is’ due 
to the fact that the linear search often finds a minimum in 
the infeasible region. However by checking for the 
feasibility of the point generated and by taking suitable 
action when it is infeasible, the method is readily adapted 
to the solution of these problems. Its performance is 
generally so superior to that of the Golden Search that it 
virtually precludes the use of that technique except under 


certain circumstances (See section 2.7.1 and 2.7.2). 


The effect of retaining the old H instead of setting it 
equal to the identity matrix at the beginning of each SUMT 
iteration is dramatic. Not only is the rate of convergence 
2 to 4 times faster but sometimes a solution is found where 
reinitializing the metric H results in no solution. 
Therefore, in general, the metric 4H should not be 
reinitialized unless during the solution of some particular 
problem the accumulative error results in H becoming either 


Singular or inaccurate. 


Rosen's GP was applied to the solution of TP13. The 
SD, PARTAN, CG and DFP implementations of the algorithm were 
applied to the solution of this problen. in addition the 
SUMT algorithm was used to solve this problem in order to 


provide a basis for comparing the two different techniques. 
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The LTesults are given in Table 3.11. 


All the implementations of the GP algorithm provided 
accurate solutions with a small execution time. For a 
relatively well-behaved objective function such as in TP13, 
CG and PARTAN performed the best with DFP possessing the 
largest execution time and largest number of iterations. 
The larger executon time is due to the more complicated 
algorithm while the much larger number of iterations is due 
to the particular path selected and is peculiar to the 
problem and not typical. The SUMT algorithm behaved very 
poorly on this problen. The PARTAN and CG algorithms 
generated solutions but they were inaccurate with more than 
14GAserror in the valine of x* ‘and «f£ (x); Only the DFP 
algorithm found an accurate solution in a reasonable amount 
of execution time. In addition the DFP algorithm had to be 
restarted at each SUMT iteration; otherwise the inaccurate 
solutions of the first three SUMT iterations contaminated 
the metric H so much that it was unable to converge to the 


true solution. 


The advantage of a direct method such as GP in solving 
constrained minimization problems over an indirect method 
such as SUMT was also apparent: GP displaying much faster 


convergence with fewer function evaluations and a smaller 


execution time. 
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3.9.5 Conclusions 


See SS aS sae 


The simple penalty and barrier methods are capable of 
producing approximate solutions and are simple to implement. 
However, some care must be exercised in choosing a 
minimization technique as some may perform well on a 


particular problem and others may not converge at all. 


While the execution time of the SUMT algorithm is much 
higher than that of the true penalty and barrier methods, it 
is much more likely to converge to the solution and the 


solution produced is much more accurate. 


In ircumstances where it can be used, a direct method 
of solving constrained problems such as gradient projection 
likely has superior performance characteristics over 
indirect methods employing penalty/barrier concepts such as 


SUMT. 
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CHAPTER LV 


4, Organizing the 


(Sa = ee ae —— mm a a a a a a SS SS ee = 


Programming Solution 


In this chapter the methods of mathematical programming 
will be applied to the optimal control of selected dynamic 
systems. In order that these optimal control problems can 
be solved, specific formulations of the problems are 
required which will be developed in the following sections. 


These will serve as the basis for further discussion. 


ee Se SS SS Se ee ae SS 


In formulating the optimal control problem we will deal 
first with continuous-time dynamic systems. From these we 
will ccnsider a discrete approximation enabling a solution 
by computer techniques. The dynamic behavior of each system 
will be described by a set of ordinary non-linear first- 
order differential equaticns of the form 

x=f (x,u,t) 
where: 

xcE” is an n-dimensional state vector; 


is an m-dimensicnal control vector; 


f(xnu pe ek. is an n-dimensional non-linear vector 
LUNCELON 3 
x= ax/de 


In addition, one will be given the boundary values for the 
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system such as the initial and/or final state. 


The control action consists of bringing the system from 
some initial state x(t.) tc some final state x(t.) in such a 
Manner that a specified performance index is minimized 
without violating the dynamic equations or any constraints. 


The general form of this performance index is 
te 
J(x,u,t)= Jreeat) dt+ P(x (te)) 
to 
In many practical applications constraints are imposed on 
the state and cecntrol vectors. These constraints are 
assumed to be cf form 
g (x,u) <0 
h (x,u) =0 
Some of the most commonly used cost functionals are: 


be 


(a) Minimum time: J= { at=t,-t, 
bo 
tf 
(b>) Minimum effort: a= | ase dt 
fo 
te 
(c) Minimum energy: a= | ue pat 


to 


4.2 Solution of the Optimum Control Problem 


The main theoretical approaches to optimal control 
synthesis are based on the Calculus of Variations [17] and 
its derivatives: Pontryagin's Maximum Principle [32], and 


Dynamic Programming [4]. 
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4.2.1 Calculus_of Variations [17] 


The classical calculus of variations, developed in the 
seventeenth century, can be used to solve the optimal 


control problen: 


te 


"Minimize J= Pexceareters | roc) sce) tat 
fo 


subject to x=f(x(t),u(t),t)" 


With 2: (fs) kand7or x{t7) (constrained. 
Define the Hamiltonian H as 
Be (Ela (cipp(e) el L(xkt) , Ute ele) Eye ete 


Then standard theory gives the Euler-Lagrange equations 


== Ot 
Ox x ox 
A(t, ) oi ue 
dH=0 
ee 


In solving these equations one has less than a full set 
cf boundary conditions specified at each end-point, 
resulting in a two-point boundary value problem (TPBVP) 


which is generally difficult te solve. 


Pontryagin's minimum principle (32] consists of a set 
of locally necessary conditions for optimality. This 
involves the use of variables called costate or adjoint- 


system variables, (A,, i=1,---,n), Similar to the Lagrange 
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multipliers. 


There exists an adjoint vector A(t) such that: 


x* (t) =0H 
y(t) 
X(t) =- 2H 
doc (€) 


TERuU*s*(L) isSocothe. optimad “control -and» x*(t)> Pheroopti nal 
trajectory then 


H (x*,u*, \*,+t) <H (x*,u, \*¥,t) 


This ccndition replaces dH/du=0 for uéU. This also results 


in a TPEVP. 


Dynamic programming is based on the observation that 
"An optimal jpolicy has the property that whatever the 
initial state and initial decision are, the remaining 
decisions must constitute an optimal policy with regard to 
the state resulting from the first decision" (Bellman's 
principle of optimality). The method of dynamic programming 
was first developed by Bellman for discrete systems, and in 


the continuous form by Hamilton, Jacobi, and Bellman. 


The major computational disadvantage of dynamic 
programming is that its recursive relations require an 
excessive amount of computer memory though this requirement 


has been partially reduced by other authors. 
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In one method of fcrmulating the optimal control 
problem as a mathematical programming problem, the 
continuous control problem is represented as a discrete 
control problem. The continuous optimal control problem is 
taken to be: 

“pind the control u*(t) U such that: 

X(t) =£ (x(t) ,u (t) ,t) 
with boundary conditions (not ae fixed) 
X{to)=X,7 X(t.) =Xe, 
follows the trajectory x*(t) X and minimizes the cost 


Function 
7. 


T(x,u,t)= D(x (t,))+ re (t) ,u(t) ,t) at" 


to 

First the time interval [t,,t,] is partitioned into N 
intervals, not necessarily all equal. The time points under 
consideration become t,,t) reed ty where ty=t2. We mast now 
approximate the state equations by difference equations. In 
doing sc, we shall use the simplest approximating difference 
equation: 

x(t, +At,) = x(t.) +£ (x (tk) -U (tk) -tk) Atk 
For brevity x(t.) is written as x(k). Then each state can 
be expressed as 

x(k+1) = x(k) +£ (x (k) -u(k) sk)AG Yk Ee (0,N-1] 


Thus the control vector u(t) becomes a sequence 


. 7 
lsh tate YBa tie, ie 

a a 
uit ,.eehdum) | pl Sumpage 
eeeT Iai -4 <6 ae 


i estiope dota & 


: tn Pivcdagy i pe) ed (og dene wae, 


ASG Seine) (rita Da ) ~~ 

y ae ate 
nage zi tor) dire RES vaubstee ari 7 
| via : Os San 
(i S® (ye) F ane i 
WAps ‘> cagseeak ier’ feite deere bese. ody sit 
| | 2 
4 7 - 


CAN et Ure d C2! rp +6 -iegnein ei : 
(erent - 


OTP) (OMNSESey «tae AORTA: eae ahs farts — 
> te 
reneyv Orit % ‘nt? A HU: aL #2 4.16064 Sad +04  tfavaerak. 
~~ a 
©. . G14 (eee é2 ‘iene eets 2 @8ouee tot setebieaeg 
. — 
(tm Meg sae TH ys naw Cue stame Me eeea nee 7 
eee Habe peg fe EA oet Gee he Poe ie 


| sberanpe 

noth he Rees Cay tye f ya be i (2s ae ~ 

a tec. Mao 1 iN hehe h te SB we tyerhx etiveepe 
el 7 - 

Pi-ayey: rer reer CO rg AR TeV at yw * Gh edpe 

Pailin ARR BOT . —a bit).  aaei* Rectat edt 

ee op bipid | 


- 


85 


{u(0),u(1),...,u(N-1)} and the state vector x is a sequence 


{x (0) -xX (1) ,-..,x(N)}, uae > 20, 


The integral cost functional can be represented as the 


limit of a summation. The system then becomes 


"Minimize 
Vv-! 
J(x,u,k)= ria Gx) x D2 eck)». 0) 1) 


Ssubjece (to. Lim 17 At. (x (k*4) == (Kk) )=f (Kk) Fa) 
At, ?o 


h (x (k) ,u (k) ) =0 u(kjyéeu 


oe cc x (k) €X 


where At, Sa es 


AWGN G x(t. y=xe 7 x (c j=" 


Thus the continuous optimal control problem is a 
mathematical programming problem of infinite dimension. 
However in practical applications the number of variables 
must of necessity be finite and this number must be small 
enough that computer implementation becomes practical. The 
above problem is then reformulated as 

Vv—| 
"Minimize J(x,u,k) = d (x (N)) + \ "1 (x (k) 0 (K) pk) Ate 
K=0 
subject to x (k+ 1) =x (k) + At, £ (x (k) ,u (k) ,k) 
g (x (k) ,u (k)) <6 
h (x(k) ,u (k)) =0 
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and x(C)=x,, x (N)iexe" 
While other more elaborate approximations of the integral 


can be used, this formulation is the simplest and most 


convenient to use. 


The problem can now be solved using mathematical programming 
techniques such as penalty methods and Rosen's gradient 


projection. 


If Rosen's gradient projection method is used and the 
state difference equations are non-linear then a technique 
employed in the method of quasilinearization must be used: 
The state equations are linearized about a nominal state- 
control history and then a sequence of linearized problems 
is Ysolved* from an ° initial ‘state-control history (ustally 
guessed) [48]. However in this study, for simplicity only 
linear systems are solved by Rosen's method and for further 
ease of programming all time intervals At are considered 
equal. The formulation given above will be applied to the 


solution of a sample problem in Chapter 5. 


Programming Considerations 


=== — 


When solving the optimal control problem as a 
Mathematical programming problem, in the format given above 
each component of the state or control vector is a separate 
variable at each discrete time point t, and in many problems 
the time points themselves are variables. This gives N time 


variables, n(N+t1) state variables, and mN control variables 
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fOrea total jofs Nidtntnj +m variables: liv Gany boundary 
conditicns such as x(0) or x(N) are given then the number of 
variables can be reduced accordingly. It should be kept in 
mind always that the above formulation is an approximation 
and the results should be interpreted appropriately. In 
some circumstances, such as high-frequency oscillatory 
systems this approximation will not work. It should be 
noted also that more elaborate difference schemes may be 


used to improve the accuracy of the sclutions. 


The epsilon method is another approach to the 
optimization of dynamic systems which avoids the explicit 


solution of the dynamic equations. 


Given the problem 
te 
"Minimize | born, tat 


subject to x=f (x,u,t) 

given x(t,)=x, and x(t,)=x-" 

the epsilon method seeks to minimize the functional 
te te 


S(xeayt se) = 26] NiR£ 000 2qdttyiai (x7, tyat 
to te 


as é€ 0. 


An iterative approach is used in which the boundary 


conditions are always exactly satisfied but where a dynamic 
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error term (the norm term) in the system dynamics exists. 

This error term is reduced with each iteration until the 

System dynamics are satisfied. The method embodies some of 

the features of both parameter optimization by differential 

approximation and parameter optimization by the penalty 
function method. The following are the three essential 
points of the epsilon method as given by Taylor and 

Constantinides [39] from Balakrishnan [2,3]: 

(1) The epsilon technique offers a non-dynamic formulation. 
The dynamic equations of the system are not explicitly 
solved; 

(2)As the epsilon parameter approaches zero, the epsilon 
method yields a statement of Pontryagin's minimum 
pranciple « Thus the above functional provides a 
sequence of trajectories and controls that can be made 
arbitrarily close to the true optimum as epsilon 
approaches zero; 

(3)Not only is ne J(x,u,t,é) likely to have a solution but 
the epsilen formulation may have a solution even when 


the solution of the original problem does not exist. 


In order to represent the epsilon functional so that 
x(t) is not identically f(x,u,t), use is made of the 
Rayleigh-Ritz method. Also, since the Rayleigh-Ritz method 
is a method of parameter optimization, the functional 
minimization Originally required becomes an ordinary 


minimization. 
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The Rayleigh-Ritz Expansion [3,39] 


SS Se Se a a a a ee a Ss 


The method as used here consists of expanding the state 
variables and control variables in terms of an infinite 


sequence of functions 
Pp. (t) gv T= es anoo 
Co 
So that x; =p, (t) +) a; p; (+) j=1,000,0 
c=! 


where p,(t) satisfies any boundary conditions and where the 


p; (t) vanish at the initial and terminal times. 


Thus the cost is a function of the variables a+. which 


4 
are determined through ordinary mathematical programming 
techniques. This technique can also reduce the total number 


of variables used in the problem. In this study the 


funetions (sin (ivt/t,)) are chosen as the p; (t). 


Other Ritz formulations can be used. These would allow 
the use of Chebyshev polynomials,cosines, etc., in place of 
the sine functions’ used. However, because of the 
modifications required in the programs and because of a lack 
of time and resources their implementation is left as a 


subject for further study. 


Balakrishnan's epsilon method is applied gee the 


solution of a test problem in Chapter 5. 
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CEAPTERV 


sie Case Studies 


The purpose of this chapter (is tot appiy’ the 
mathematical programming techniques discussed in chapters 2, 
3, and 4 to the soluticn of two selected control problems. 
The examples presented were chosen to demonstrate these 
techniques and as such are relatively simple and straight 
forward. Each problem is solved by more than one approach 


in order to illustrate and compare methods. 


The first sample control problem is a linear system 
With a quadratic cost function and initial and final states 
specified. It will be referred to as the minimum energy 


regulator (MER) problem. It consists of: 


te 
"Minimize J= { u2 (t) dt 
E 


(o) 
subject to Grey 
=u 


with the boundary conditions 


t,=0.9 t, =2.20 


plus the additional constraints 
g (x,t) =x, +3x,+0.02t =1o; 0 


-1 <u <1 
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The analytical solution for this problem without g(x) can be 
easily found [19], but otherwise a computer solution is 
required. 
In order to. solve this problem using mathematical 
programming techniques it is discretized in the following 
Manner. Divide [t,,t,] inte m equal intervals of length At. 
The problem is then reformulated as: 
m-' 
"Minimize J(x,u,t) =At) u2(4) 
i=o5 
subject to x Gtly=x 0G) +Arx (a) t=O; 1 00a O| 
x(t) =x.) Aen) 1=0 7 Typo 
g(x (ip t (2))) Hun P42) 432s (a) 402028 (1) 150 
Ss (a) as) 
where At=t, /m 
£qayaide® 


Pei cues jie (2) 


The above MP problem now has 4m+3 variables and 5mt+t3 
constraints. For large m this would be expensive to solve. 
Using the difference equations and boundary conditions we 
Can solve for the x(i) explicitly in terms of the wu(i). 
This results in: 

i=! 
x, (4)=x, (0) +At) x, (J) A= lip oerey 
Hee 
i-l 


70(2) =x, (0) t4¢/ 0) (a) f= 4 erery i 
' oO 


| = 


ru (2) =x (0) +idt 
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Substituting for the X, (Jj) in the equation for x, (i) yields 


X, (1)=x, (0) + Atx, (0) 


jet 


i- 
x, (1) =x, (0) +idtx ,(0) +dt2) ,) a (ky 1=2 7 suis 


j=l kso 


t-! 
=x, (0) +iAtx, (0) +At2 ) \(4-k-1) u(k) i 


K=O 


1 gece og m 
SUDStIcUting into g(x ())<0° yields 
i-| Tea 
spn ig gag AN ae (0) +At2)> \(4-k-1) 0 (k) #3 (x, (0) + At) u (Kk) ) 
K=0 k=0 
#0.02(x,(0)+iAt)-1.5 <0 LEW ete pt 


Regrouping terms and substituting the BC's, then, 
i-| 
g (x (4)) =) (At? (i-k-1) #3At Ju (k) #0.02 (44t)-2.5 <0 


k=0 
P= lee pil 
The equalities x, (m)=0 and x,(m)=0 reduce to -x, (m)<0 and 
x, (m) <0 as the equalities will hold if the energy is 
minimized. The inequalities 
=k; (m) 50 and: =x, (m)<0 


can be expressed as 


M-x 


-At2) (m-k-1) u(k) +1 <0 
k=0 
m— | 
and At > ‘a (k) <0 
k=O 


Thus the problem has been reduced to an MP problem of am 


Variables with 3m+2 constraints. The reduced probiem is 
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restated briefly: 


m-| 
"Minimize J=At) ‘u(iy 2 


ize 
subject to u(i)-1<0 DO i eree gin 


i] 


[At2 (i-k-1) +3At Ju (k) 0.02 (At) -2.5 <0 sf 


) 


Zi paterere al 


KR 
1 


Mez 
- dt2 )|(m-k-1) u(k) +1 <0 
k=0 


ml 
At x (k) <o" 
k=0 

This problem is now in a form that will allow it to be 
Solved directly by pregrams implementing methods such as 


Kosen's Gradient Projection method or Fiacco and McCormick's 


SUMT. 


Since the program does not contain a procedure for 
generating a feasible point from a non-feasible one and 
because of the difficulty in choosing a feasible point for 
non-trivial m, a slack variable, S, was introduced into the 
constraints: 

Ep u(t )j=0.5 =O pipes ap ea a | 
u (2,)-=-0..5 B= /2ij eer a 


then all the constraints are satisfied except 
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m-2 
-At2>!(m-k-1) u(k) +1 <0 
k=0 


if the Slack variable S, with S°=1, is added to at, i 


becomes 
m-2 
- (\t2 eat u(k)-S+1 <0 
K=0 


and we have an initial feasible point. In order to have an 
optimal feasible soluticn tc the original problem, the value 
of S must be 0 at the constrained optimum. This can be done 
by modifying the cost functional J by adding the penalty 


term 500S2. The cost functional thus becomes 


m-| 
J=At Du (4) 2450082 
=o) 


The above method for obtaining an initial feasible 
point was adopted except that, since logarithmic barrier 
functions were used, the two constraints incorporating the 


slack variable were modified to prevent negative arguments: 


m—-Z 
- Ate) \(m-k-1) u(k) -S2+1 <0 
k=O 


m ~ 


- At > 'u (k) -S2<0 
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Figure 5.2 MER optimal trajectory 
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5.2 Numerical Resuits 


The results of applying Rosen's GP and SUMT to the 
solution of the MER problem are presented in Table 5.1. The 


optimal trajectory and control are given in Figures 5.1 and 


5 e 2 e 

fae he eo , Wan ee a ee ae a ee 
{TEST|PROB| SD | CPTAN{ CG | DFP | 
fon nf 
{GRAD| IT | 30 { 25 | 2h | 2061 
{PROJ| #F | 39 | 2s | oh | oral 
| (TIME116.5894 14.822414.604420.231| 
{- +} ---- {=== --4------ }------ }------4 
| , it | # | 138*{ >280*| 159 | 
([SUMT[ #F | 1 848 1>1296 | 634 | 
{TIME} {58.8081>120.0158.759| 
(ee? ee ee ae wl ea es | 


m=30 time steps, r=1.0, c=5.0 
# - NO SOLUTION ATTEMETED 
* - INACCURATE ANSWER 


Table 5.1 Numerical results: MER 
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omments 


All four implementaticns of GP (SD, PARTAN, CG, and 
DFP) solved the problem quickly and accurately for m=30 time 
intervals. Conversely the SUMT algorithm had a great deal 
more trouble, with CG being unable to produce a solution in 
a reascnable amount of execution time and the solution 
provided by PARTAN being highly inaccurate. Only DFP 
provided a reasonably accurate SUMT solution in a practical 


length cf time. 


As expected for the GP solution, SD had the largest 
number cf iterations and DFP possessed the largest execution 
time, the objective function being a well-behaved quadratic. 
In general, the behavior of the accelerators was as detailed 
in 3.9.4 except that DFP had the smallest number of 
iteraticns. Again the advantage of using a direct method as 


compared to an indirect method was apparent. 


The solution provided by the DFP SUMT algorithm 
required 634 function evaluations and 3 times the amount of 


executicn time. 


Additionally it was much easier to calculate 
analytically the function and and gradient required for the 
GP algorithm than it was for the SUMT version (another 


common advantage of this method of solution). 


Conclusion 
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These results suggest that this implementation of GP 
incorporating several accelerators is an effective tool for 
the mathematical programming solution of control problems 
with linear constraints and in general is superior for this 
type of problem to simple implementations of penalty-barrier 


algorithms. 


5.3 Earth-Mars Orbit Transfer Control Problem 


———_—— me ee a ee eS ee ee eee ee ee ee 


The second sample control problem is a non-linear 
system with a minimum time cost functional, which will be 
solved using the epsilon technigue [39]. This problem is 
the normalized Earth-Mars orbit transfer solved previously 
in a paper by Moyer and Pinkham [28] and by Taylor and 
Constantinides using the epsilon technique. The orbit 
transfer is (to be performed using minimum fuel, but for a 
constant-thrust ion jet this is equivalent to minimum tine. 
The control variable is the thrust angle 9. The system 
equations are: 


oe eee 
X, =x2-( M/x2)+(Tsin6)/(m, +mt) 


x, =~ (%,X, /X, )+(Tcos 6) / (m, +t) 
The boundary conditions and constants are 
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The epsilon cost functional is 
te 


Ee 
scorer= 11/28 [as-eoe.8.e) uate | 
to b 


dt 
The geometry of the problem is given in Figure 5.3. 

The state variables and control variable will be expanded by 
the Rayleigh-Ritz (RR) method. The easiest and most obvious 


function to use in the expansion zs the function 


sin(itt/t,). The discrete representation is now: 
M 
Reh tie e hi ee 1) yay Sin(; ) 


d=! 
M 


<= (ke = Ke N-1) At + ji a costa.) 
in = (Kig Xv 70-1} 2 Pi ( 0; 


™M 
O,= 8+ ( 8,- 8) (n-1)7 (N-1) + > \b, sin ();) 
i 
where 
N=total number of time points 
n=time point under consideration 


M=number of terms in the Rayleigh-Ritz expansion 


(= Jr(n-1)/ (8-1) J 


ea (8-1) Ae 
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Figure 5.3 Geometry of Earth-Mars orbit transfer 


The parameters ar 1dr O84 Oy and At are the variables to be 


solved for. 


Let the epsilon cost functional he expressed in the 


form 


to 6)\= 61S 


where 6 is a (3N+#1) column vector given by 


S(z)=E [At Ze Sin Ve AVAE tend, VOX 8h, 0, (OH 
1) A954 


and 
OX in Sx ~Xon 
=Koq — (K2, /Xm )+/Yx2 )- (Tsin O,)/ (m, +m (n-1) At) 


(mar Hee ee) 7k Cos On, (mem (a1 Ae) 


and zis the (4M+3) column vector of variables 
z=[a,, peee gay 7ar geeeeaan 9 a2, peooe garay roy Parry ee Or Oc, At J 
(The above formulation, with corrections, is from the paper 


by Taylor and Constantinides [39)]) 


The resulting problem can be solved in three different 
ways: 
(a) The Gradient Technique: This is a standard method for 
solving control problems. However Taylor and others 
[3,41,42] have been unable to obtain satisfactory results 
using this method. A solution by this method was not 


attempted. 
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(b) The Modified Newton-Raphson Method (MNR): This method is 
mentioned by Balakrishnan [3] and developed by Taylor et al 
{41] and was used by Tayler and Constantinides [39] to solve 
this preblem. A development of it will be given below. 

(c) Other Mathematical Programming Technigues: This probien 
can be solved using algorithms such as DFP, conjugate 
gradients etc. A short development of these solution 


techniques will be given below. 


Using this method “the cost functional “7(  , ) is 
expressed as 6'§ . The vector 4 is linearized in terms of 


the change of the problem variables, Az. Thus 
6=8,+F, (5) Az 


' 
g= 6,5 +25, Azedz FF, Az 


where z was defined earlier. 


F,(¢6) is the Jacobian matrix of 4. The minimization 
is accomplished by iteratively solving 
Ze 


=Z tH AZ: 


Setting the gradient of J(&,€) to zero 


en ce ae 
JpoesOae on eet? Az) FEY 


and solving for the difference Az. one gets 


Az, =-(F, (5,)" F, (5,) F1F, (8), 
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523.2 Sclution by the MNR Method and Golden Search 


A small variation of the MNR was attempted whereby «the 
stepsize was optimized for the minimum value of J. This 


Optimization was done using the Golden Search. 


a a a a a a ae se a es a ss ee ae SS ee ee ee ae iw Se a as we = 


Again the cost functicnal J( 9,€) is expressed as 65°S. 
One can use first order gradient techniques to minimize 


J(6,€) after calculating the gradient as 
3(6,€)= 55 


SALE OE 
od(6,€)= 5 06 + (38) 


The problem is now in a form readily solved by such 


techniques as DFP, conjugate gradients, etc. 
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There is a second method due to Taylor and 
Constantinides [40] of representing the problem. This 
method, with fewer variables, consists of expanding only the 
state variables by the Rayleigh-Ritz method and then 
performing an ordinary minimization for the control 
variables at each discrete point along the state trajectory. 


The problem thus reduces to 
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with the state variables represented as before. The optimal 
control is then obtained at each point along the discrete 
trajectcry by minimizing the dynamic error with respect to 


the control variables. Thus setting 


d (ix-£ (x,8,t) | 2) =0 


—— 


00 


yielding O=tan-1{x,- (x2/x, )+(/Yx2)}/{%,+ (%,X5/X, )} 


at each point along the trajectory. 


5.4 Numerical Results 


This problem was solved using the four techniques 
discussed above. The results of applying these techniques 
to the sclution of the Earth-Mars orbit transfer problem are 
presented in Table 5.2. The optimal trajectories and 
controls produced by each method are given in Figures 5.4 to 
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Figure 5.5 State trajectory X, 
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Figure 5.7 Optimal control angle 
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Figure 5.8 Cost versus iterations - NR minimization 
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Figure 5.9 Cost versus iterations - DFP minimization 
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Figure 5.10 Error versus iterations - NR minimization 
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ragube 5.11) Error versus iterations —- DEP ninimi zation 


Taylor and Constantinides! results for N=30, m=7, and 
€=10-* were reproduced for the MNR method. However in order 
to achieve similar results, the technique (mentioned by 
Taylor etal. [42]) of increasing the diagonal elements of 
the matrix [FF] by a small percentage had to be adopted. 
Without this modification the solutions did not converge and 
oscillated wildly due to the near singularity of this 
Matrix. The value of the increment was not critical in this 
case and 5% was chosen. The number of iterations was set at 
20 and the initial control(9,6,) was set to (0,217) since 
this gave reasonable results for the DFP solutions and 
computing cost was a LacLor which inhibited wide 
experimentation. The DFP algorithm was used to provide all 
other solutions presented. The other aigorithms performed 
worse, with CG taking twice the execution time to provide a 
solution and with PARTAN taking even longer and not 
converging to an acceptable solution. SSVM was comparable, 


especially if restarted. 


For the formulation with the control variable expanded 
as a kR series the MNR and DFP solutions were almost 
identical. However the DFP solution took less than half the 
execution time and 50% less memory. This is due to the 


(4M+3) x (4m+3) /2 matrix inversion required for the NR 
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solution. This situation becomes even greater in the DFP 
solutions' favour as the number of variables is increased 
because most of the cost of the MNR method is in this 
inversion. This is expensive even though advantage is taken 


of the fact that the matrix is Symmetric positive-definite. 


The method employing the Golden Search as a Single 
variable search worked and produced solutions almost 
identical (within .01%) to the DFP solutions for the (4m+3) 
variable problem. The number of iterations was half that of 
the MNR method and required only approximately 60% of the 
execution time. This accelerator technique did not work if 
the elements of the diagenal were modified as in the 
Standard MNR method, the solution times being the same in 
both cases. The tolerance of convergence chosen for the 


Golden Search was a crude 0.4. 


The Golden Search did not work for the VMNR method, 
cauSing the problem to converge to a non-optimal solution. 
The technique of incrementing the diagonal elements also was 
unsuccessful for this method, again causing convergence to 
the same non-optimal solution.’ However the "convergence" 


was more rapid. 


The VMNR method of sclution was not found to be very 
Satisfactory. The control was very sensitive to changes in 
t,, m, and € and varied significantly as these were changed. 


Also while the cost per iteration was less, a larger number 
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of iterations was required for a solution of comparable 


accuracy. 


In addition, this formulation possessed at least two 
other local minima to which the method often converged. 
This fact plus the lack of information in the article [40] 
made it impossible tc achieve comparable results. AS a 
consequence some of the conclusions reached here differ from 
those given by Taylor and Constantinides for this method 
{490]. These new conclusions are: 

(1) The convergence of the algorithn requires more 
iterations for a comparable control and accuracy, though the 
cost per iteration is less. 

(2) The dynamic error is greater resulting in a less 
accurate control, though this error can be reduced by more 
iterations. 

(2) This method created at least two other non-optimal 


solutions to which the method had a tendency to converge. 


The DFP solution of the formulation involving (smt1) 
variables was even more sensitive, the solution time for a 
comparakle cost being approximately the same as for the VMNR 
method. The VMNR solutions for this problem tended toward 
the control trajectory shown in Figure 5.7. This trajectory 
is similar to that produced by Moyer and Pinkam [28] for 
their solutions to the Earth-Mars orbit transfer problem 


produced by the Gradient and Second Variation techniques. 
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The main disadvantage of the DFP solution of these 
problems was that the soluticn time was sensitive to the 
initial boundary conditicns. The solutions times of the MNR 
method were relativehy unaffected by changes in €, initial 
total time, and initial valves of the control; on the other 
hand, the DFP sclution time was affected by changes in these 
values and solutions both much better and much worse than 
those given were achieved. The initial conditions chosen 
represented the mean values and also were a logical initial 


guess. 


The choice of the parameter had a great effect on the 
rate of convergence of the DFP solutions. For large the 
DFP algorithm converged much more rapidly than for a smaller 
walue of €. However, as in the NR case, the error in the 
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CHAPTERGVE 


6. Summary and Conclusion 


———— ee te ee ee ee ee ee ee ee ee 


This study has been concerned with the review, 
implementation and evaluation of Standard nonlinear 
mathematical programming techniques and their application to 
the soluticn of certain mathematical programming and optimal 


control problems. 


6.1 Summary of Results 


in Chapter 2 the metheds used in solving unconstrained 
problems were reviewed and implemented practically. These 
computer implementations were applied to the solution of a 
family of test problems. The solutions were then compared 
and the relative merits cf the algorithms were discussed. 
The conclusions drawn agreed generally with those of other 
authors in this field except that. the results were often 
found to be highly dependent on the mode of implementation 
and to a lesser extent dependent on the problem under 
consideration. From this it can be inferred that there is 
no way cf determining a pricri which method will be best for 
solving a particular problem, but only that some methods are 


more likely to prove more successful than others. 


The solution of constrained problems was the topic of 
study of Chapter 3. In this Chapter the solution by penalty 
and barrier methods, including SUMT, was discussed and then 


applied to the solution of a set of test problems. The 
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solutions of these problems by the different methods using 
the minimization techniques of Chapter 2 were then compared. 
This comparison revealed the defects inherent in these 
methods and the weaknesses of first order gradient 
minimization techniques when applied to the soluticn of 


penalty constrained problems by these techniques. 


Rosen's gradient projection algorithm was described 
. along with an implementation of accelerating techniques 
which improved its rate cf ccnvergence. These accelerators 
consisted of the PARTAN, conjugate gradient and DFP 
algorithms. The accelerators, particularly CG and PARTAN, 
were found to improve substantially the rate of convergence 
in comparison with the original steepest descents direction 


generation. 


The performance of gradient projection was also 
compared to that of SUMT and it was shown that for a problem 


with linear constraints GP is more efficient. 


Chapter 4 reviewed two methods of transforming the 
optimal control problem into a mathematical programming 
problem. The first method consisted of discretizing the 
control problem and approximating the state eguations by 
difference equations. The second method discussed was 
Balakrishnan's epsilon technique, a method of optimization 


which avoids the explicit solution of the dynamic equations. 


Chapter 5 consisted of applying the mathematical 
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programming techniques of Chapters 2 and 3 and the methods 
of Chapter 4 to the solution of two selected control 
problems. These problems were the minimum energy regulator 
problem and the Earth-Mars orbit transfer problem. MER was 
solved ty means of Rosen's gradient projection and SUMT with 
the twe methods being ccmpared. The Earth-Mars orbit 
transfer problem was solved by means of a modified Newton- 
Raphson technique and by the DFP algorithm. The two methods 
of solution were compared with the superiority of the DFP 
algorithm for this type of problem (possessing a large 


number of variables) being demcnstrated. 
6.2 Conclusions 


In the course of this study, possible areas of further 
investigation became apparent. Relatively little work has 
been done in applying mathematical programming to the design 
of optimal control systems. More particularly, mathematical 
programming techniques could be used in applying 
Balakrishnan's epsilon technigue to a wider range of control 
problems, especially those whose cost function is something 
other than minimum time. More work could also be done on 
extending the types of Ritz expansions used in the epsilon 


technique, especially to the use of expansions other than 


the sine series. 


Ancther area of interest is on-line use of the epsilon 


technique in the generaticn of controls for dynamic 


’ tt Ldow 


ee ed Coe Le ced | oat can her 
‘ Ps rt epee & he ' BD 04 ern 2 


im . : Ty ‘ergtyd's hing ; ' a Bhodten 


f etre { Vi em, | banana: eae! 
2 a ee 


y) 
MacimUm ps ety Hennuds aren! apy 


J 
t) he T ss ? aS ig 


ayn wud aobrepiravad 


Logitéeedtaa Oct @higyauil ag dead 
; , 


‘ By } np ue ei eiwy is 
\ Lv LT 4) "6° 
A) Ao i Tie ee a 17 = 
Pie! jgisds aor age, Gucd sy hl adie jw 
pi iaay 3 34.) BL S- -eehsy itll =e 


bee ituberuw 2eid eo & 


re © ee 


12a 


processes. This would entail studying the accuracy of the 
controls generated, the speed of solution and the minimum 


memory requirements necessary for acceptable performance. 
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